// do /Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Code/plotting/plot_Fox_contour_3dim_bs.do
//Options: discrete 44000 16 1

pause on


global  pathresults = "/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_estimation"
global  Est_pathresults = "/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_estimation/demand_est_grID"

local FOX_model = `"`1'"'   // discrete, 
//Robustness: discretev5, discretev6, discretev7, discretev8, discrete_3tau 
local subsample = `"`2'"'   // Use 44000 
local inc       = `"`3'"'   // 16 1  2 3 4 5 6
local set_seed  = `"`4'"'   // 1  10 3 6 4 3 7 




insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear


 forvalues i=1(1)16 { 
	capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(theta eta)
 
//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 eta_avg theta_avg, i(theta eta) j(bs)

keep if weight_100>0
drop if weight_100==. 
sort eta theta
by eta theta: egen mean_weight=mean(weight_100)
by eta theta: egen sd_weight=sd(weight_100)
by eta theta: egen count_weight=count(weight_100)
gen  se_weight=  sd_weight/(count_weight^0.5)
gen ICl_weight=  mean_weight - 1.96* se_weight
gen ICh_weight=  mean_weight + 1.96* se_weight

preserve 
collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(eta theta)

	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order eta theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort eta theta
outsheet using $pathresults/plot_jointpdf_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 

save $pathresults/plot_jointpdf_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs, replace 



restore

collapse(firstnm) eta_avg theta_avg,by(bs)
	gen eta_avg_sd = eta_avg
	gen theta_avg_sd = theta_avg
	gen eta_avg_n = eta_avg
	gen theta_avg_n = theta_avg
collapse(mean) eta_avg theta_avg (sd) eta_avg_sd theta_avg_sd  (count)	eta_avg_n  theta_avg_n
	gen eta_avg_se = eta_avg_sd / (eta_avg_n^0.5)
	gen theta_avg_se = theta_avg_sd / (theta_avg_n^0.5)
	keep eta_avg theta_avg eta_avg_se theta_avg_se
outsheet using $pathresults/means_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 


//theta
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(theta)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta) j(bs)

keep if weight_100!=0
sort theta
by theta: egen mean_weight_theta=mean(weight_100)
by theta: egen sd_weight_theta=sd(weight_100)
by theta: egen count_weight_theta=count(weight_100)
gen  se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 


//Process raw FKRB results to plot the marginal eta:
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(eta)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

drop eta_avg_tmp*
//drop theta_avg_tmp*

reshape long weight weight_100 eta_avg, i(eta) j(bs)

keep if weight_100!=0
sort eta
by eta: egen mean_weight_eta=mean(weight_100)
by eta: egen sd_weight_eta=sd(weight_100)
by eta: egen count_weight_eta=count(weight_100)
gen  se_weight_eta=sd_weight_eta/(count_weight_eta^0.5)
gen ICl_weight_eta=  mean_weight_eta-1.96* se_weight_eta
gen ICh_weight_eta=  mean_weight_eta+1.96* se_weight_eta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(eta)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order eta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort eta
outsheet using $pathresults/plot_marginaleta_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 



//Process raw FKRB results to plot the marginal theta:
//theta: [<0, 0-0.75, 0.75-1.25, >1.25]
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen theta_4q =. 
 replace theta_4q = 1 if theta < 0
 replace theta_4q = 2 if theta >= 0 & theta<0.75
 replace theta_4q = 3 if theta >= 0.75 & theta < 1.25
 replace theta_4q = 4 if theta >= 1.25

 collapse(sum) weight*, by(theta_4q)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta_4q) j(bs)

keep if weight_100!=0
sort theta_4q
by theta_4q: egen mean_weight_theta=mean(weight_100)
by theta_4q: egen sd_weight_theta=sd(weight_100)
by theta_4q: egen count_weight_theta=count(weight_100)
gen se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta_4q)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_4q_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 



//Results reported in the main table where cut-offs for theta are based on the following discount rate: below 2% 2-12% 12%-very large, wrong sign
//theta: [<0, 0-0.75, 0.75-1.25, >1.25]
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen theta_4q =. 
 replace theta_4q = 1 if theta < 0
 replace theta_4q = 2 if theta >= 0 & theta < 0.6202
 replace theta_4q = 3 if theta >= 0.6202 & theta < 1.2825
 replace theta_4q = 4 if theta >= 1.2825

 collapse(sum) weight*, by(theta_4q)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta_4q) j(bs)

keep if weight_100!=0
sort theta_4q
by theta_4q: egen mean_weight_theta=mean(weight_100)
by theta_4q: egen sd_weight_theta=sd(weight_100)
by theta_4q: egen count_weight_theta=count(weight_100)
gen se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta_4q)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_4q_2_12_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 


pause


//Compute E[m] using heterogeneous r.
//
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen rho_5=1/(1+0.05)
 gen LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5)
  
 gen rho_2=1/(1+0.02)
 gen LLC_2_18=rho_2*(1-rho_2^18)/(1-rho_2)
  
 gen rho_12=1/(1+0.12)
 gen LLC_12_18=rho_12*(1-rho_12^18)/(1-rho_12)
 
 gen m_2=LLC_2_18/LLC_5_18
 gen m_12=LLC_12_18/LLC_5_18
 gen theta_m_adj = 1
 gen theta_m_2_v    = theta*LLC_5_18/LLC_2_18
 gen theta_m_12_v   = theta*LLC_5_18/LLC_12_18
        
 replace theta_m_adj=theta_m_2_v if theta > m_2
 replace theta_m_adj=theta_m_12_v if theta < m_12		
       

 collapse(sum) weight*, by(theta theta_m_adj)
 
//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
 
 capture gen theta_m_adj_avg_tmp`i' = theta_m_adj*weight`i' 
 capture egen theta_m_adj_avg`i' = sum(theta_m_adj_avg_tmp`i')
 
}  

drop theta_avg_tmp*
drop theta_m_adj_avg_tmp*
drop weight*

reshape long theta_avg theta_m_adj_avg, i(theta) j(bs)
	   
collapse(firstnm) theta_avg theta_m_adj_avg,by(bs)
	gen theta_avg_sd = theta_avg
	gen theta_m_adj_avg_sd = theta_m_adj_avg
	gen theta_avg_n = theta_avg
	gen theta_m_adj_avg_n = theta_m_adj_avg
collapse(mean) theta_avg theta_m_adj_avg (sd) theta_avg_sd theta_m_adj_avg_sd  (count)	theta_avg_n  theta_m_adj_avg_n
	gen theta_avg_se = theta_avg_sd / (theta_avg_n^0.5)
	gen theta_m_adj_avg_se = theta_m_adj_avg_sd / (theta_m_adj_avg_n^0.5)
	keep theta_avg theta_m_adj_avg theta_avg_se theta_m_adj_avg_se	   
	   
outsheet using $pathresults/mean_theta_theta_m_adj_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 
